#include "psoml_class.h"

void Psoml::SetParametersFromFile(char *path)
{
	ifstream argumentFileStream(path);

	if(!argumentFileStream)
		throw new exception();

	while(!argumentFileStream.eof())
	{
		string line;
		argumentFileStream >> line;

		string argument = line.substr(0, line.find("="));
		string value = line.substr(line.find("=") + 1, line.length());

		if(argument == "--processorcount"){
			processorCount = atoi(value.c_str());
			swarmSize = processorCount * 128;
		}
		else if(argument == "--outputfile")
			outputFileStream = new ofstream(value.c_str());
		else if(argument == "--maxgenerations")
			generationMax = atoi(value.c_str());
		else if(argument == "--currentcoefficient")
			swarmConstants.w = atof(value.c_str());
		else if(argument == "--neighborhoodsize")
			neighbourhoodSize = atoi(value.c_str());
		else if(argument == "--globalcoefficient")
			swarmConstants.z = atof(value.c_str());
		else if(argument == "--neighborhoodcoefficient")
			swarmConstants.y = atof(value.c_str());
		else if(argument == "--localbestcoefficient")
			swarmConstants.x = atof(value.c_str());
		else if(argument == "--lengthdeltamax")
			lengthDeltaMax = atof(value.c_str());
		else if(argument == "--lengthdeltamin")
			lengthDeltaMin = atof(value.c_str());
		else if(argument == "--lengthmax")
			LengthMax = atof(value.c_str());
		else if(argument == "--lengthmin")
			LengthMin = atof(value.c_str());
		else if(argument == "--randomvectormax")
			randomVectorMax = atof(value.c_str());
		else if(argument == "--randomvectormin")
			randomVectorMin = atof(value.c_str());
		else if(argument == "--topology")
			phylogeny.CreateTreeFromFile(&value);
		else if(argument == "--sequencelength")
			phylogeny.sequenceLength = atoi(value.c_str());
		else if(argument == "--sequence"){
			if(value.substr(value.find('|') + 1, value.length()).length() == phylogeny.sequenceLength)
				phylogeny.sequences[value.substr(0, value.find('|'))] = value.substr(value.find('|') + 1, value.length());
		}
		else if(argument == "--stationaryprobability"){
			switch(value.substr(0,value.find('|'))[0]){
				case 'G':
					phylogeny.stationaryProbabilities[0] = atof(value.substr(value.find('|') + 1, value.length()).c_str());
					break;
				case 'A':
					phylogeny.stationaryProbabilities[1] = atof(value.substr(value.find('|') + 1, value.length()).c_str());
					break;
				case 'T':
					phylogeny.stationaryProbabilities[2] = atof(value.substr(value.find('|') + 1, value.length()).c_str());
					break;
				case 'C':
					phylogeny.stationaryProbabilities[3] = atof(value.substr(value.find('|') + 1, value.length()).c_str());
					break;
			}
		}
	}
}

void Psoml::Run()
{
	vector<float> times;
	//seed random number generator with time
	srand(time(NULL));
	//create default kernel launch sizes
	dim3 singleThread(1,1,1);
	dim3 blockSize(128,1,1);
	dim3 gridSize(processorCount,1,1);
	//create device/host variables
	cNode* d_cNodeArray;
	cudamalloc(&d_cNodeArray, sizeof(cNode) * phylogeny.length);
	pNode* d_pNodeArray;
	cudamalloc(&d_pNodeArray, sizeof(pNode) * phylogeny.length * swarmSize);
	int* d_globalMax;
	cudamalloc(&d_globalMax, sizeof(int));
	
	float *d_likelihood, *h_likelihood;
	cudamalloc(&d_likelihood, sizeof(float) * swarmSize);
	h_likelihood = new float[swarmSize];
	memset(h_likelihood, 0, sizeof(float) * swarmSize);
	cudamemcpydevice(d_likelihood, h_likelihood, sizeof(float) * swarmSize);

	float* d_likelihoodBest;
	cudamalloc(&d_likelihoodBest, sizeof(float) * swarmSize);
	//create and initialize neighbourhood
	int *h_neighbourhood, *d_neighbourhood;
	h_neighbourhood = new int[swarmSize * neighbourhoodSize];
	cudamalloc(&d_neighbourhood, sizeof(int) * swarmSize * neighbourhoodSize);
	for(int p = 0; p < swarmSize; p++)
		for(int n = 0; n < neighbourhoodSize; n++)
			h_neighbourhood[p * neighbourhoodSize + n] = rand() % neighbourhoodSize;
	cudamemcpydevice(d_neighbourhood, h_neighbourhood, sizeof(int) * swarmSize * neighbourhoodSize);
	//create and initialize neighbourhood bests
	int* d_neighbourhoodBest;
	cudamalloc(&d_neighbourhoodBest, sizeof(int) * swarmSize);
	//create and initialize particles - random vectors between ranges
	//seed random number generator
	pNode* pNodeArray = new pNode[phylogeny.length * swarmSize];
	for(int p = 0; p < phylogeny.length * swarmSize; p++){
		float range = LengthMax - LengthMin;
		pNodeArray[p].length = LengthMin + range * 0.001 * (rand() % 1000 + 1);
		range = lengthDeltaMax - lengthDeltaMin;
		pNodeArray[p].delta = lengthDeltaMin + range * 0.001 * (rand() % 1000 + 1);
		pNodeArray[p].likelihood = 0.0f;
	}
	cudamemcpydevice(d_pNodeArray, pNodeArray, sizeof(pNode) * phylogeny.length * swarmSize);
	
	//create timer
	unsigned int timer;
	cutCreateTimer(&timer);
	//iterate over each generation
	for(int g = 0; g < generationMax; g++){
		cutStartTimer(timer);
		//iterate over entire sequence
		for(int s = 0; s < phylogeny.sequenceLength; s++){
			//create a topology consistent with the desired site
			cNode* cNodeArray = phylogeny.cudaNodeArray(s);
			//copy topology across device
			cudamemcpydevice(d_cNodeArray, cNodeArray, sizeof(cNode) * phylogeny.length);
			//calculate likelihood
			EvaluateParticles(blockSize, gridSize, d_cNodeArray, phylogeny.length, d_pNodeArray, phylogeny.stationaryProbabilities, d_likelihood);
		}
		////update local and neighbourhood bests
		NeighbourMax(blockSize, gridSize, d_likelihood, d_likelihoodBest, d_neighbourhood, neighbourhoodSize, d_neighbourhoodBest, d_pNodeArray, phylogeny.length);
		////determine global max			
		Maximum(singleThread, singleThread, d_pNodeArray, d_globalMax);
		////update lengths and deltas
		float range = randomVectorMax - randomVectorMin;
		float3 randomVectors;
		randomVectors.x = randomVectorMin + range * 0.001f * (rand() % 1000 + 1);
		randomVectors.y = randomVectorMin + range * 0.001f * (rand() % 1000 + 1);
		randomVectors.z = randomVectorMin + range * 0.001f * (rand() % 1000 + 1);
		UpdateLengthDeltas(blockSize, gridSize, d_pNodeArray, phylogeny.length, d_globalMax, d_neighbourhoodBest, swarmConstants, randomVectors);
		times.push_back(cutGetTimerValue(timer));
		cutResetTimer(timer);
	}
	cutDeleteTimer(timer);

	averageTime = 0;
	for(int i = 0; i < times.size(); i++)
		averageTime += times[i];

	averageTime /= times.size();

	//FREE MEMORY!!!
	cudaFree(d_cNodeArray);
	cudaFree(d_globalMax);
	cudaFree(d_likelihood);
	cudaFree(d_likelihoodBest);
	cudaFree(d_neighbourhoodBest);
	cudaFree(d_pNodeArray);

	free(h_likelihood);
	free(h_neighbourhood);
}
